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Presented  by 


Abstract 

We  present  a  certified  reduced  basis  (RB)  method  for  the  heat  equation  and  wave  equation.  The  critical  ingredients 
are  certified  RB  approximation  of  the  Laplace  transform;  the  inverse  Laplace  transform  to  develop  the  time-domain 
RB  output  approximation  and  rigorous  error  bound;  a  (Butterworth)  filter  in  time  to  effect  the  necessary  “modal” 
truncation;  RB  eigenfunction  decomposition  and  contour  integration  for  Offline-Online  decomposition.  We  present 
numerical  results  to  demonstrate  the  accuracy  and  efficiency  of  the  approach. 
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Resume 

Une  methode  de  bases  reduites  certifiee  utilisant  la  transformee  de  Laplace  ;  Application  l’equation 
de  la  chaleur  et  l’equation  des  ondes 

Nous  introduisons  une  methode  de  bases  reduites  certifiee  pour  L  equation  de  la  chaleur  et  l’equation  des  ondes 
utilisant  la  transformee  de  Laplace.  Les  ingredients  essentiels  sont  les  suivants  :  une  approximation  par  bases 
reduites  certifiee  de  la  transformee  de  Laplace,  une  transformee  de  Laplace  inverse  pour  1’ approximation  de 
l’output  par  bases  reduites  en  temps  et  l’etablissement  de  bornes  d’erreur  correspondantes  rigoureuses,  un  filtre 
en  temps  (de  Butterworth)  pour  mettre  en  place  la  troncation  “modale”  necessaire,  une  decomposition  en  fonctions 
propres  par  bases  reduites  et  une  integrate  de  contour  pour  la  decomposition  Offline-Online.  Nous  presentons  des 
resultats  numeriques  qui  demontrent  la  precision  et  l’efficacite  de  l’approche. 
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Nous  considerons  une  equation  de  la  chaleur  et  une  equation  des  ondes  parametree  (par  un  parametre  /i) 
avec  une  forme  sesquilineaire  m  (de  masse)  et  une  forme  sesquilineaire  a  (de  raideur)  toutes  deux  affines 
en  le  parametre,  voir  (1).  Une  projection  de  Galerkin  est  supposee  fournir  la  solution  “de  reference” 
G  X*  par  approximation  elements  finis.  L’output  est  donne  par  une  fonctionnelle  filtre  de 
Butterworth  du  champ,  voir  (2).  Nous  reformulons  le  probleme  en  terme  de  transformee  de  Laplace  : 
pour  une  frequence  donnee  u>,  la  fonction  (u>;  n)  satisfait  A{u^  {w,  n),  v;  u>;  fi)  =  g(iui)f(v),\/v  G  . 
Dans  cette  relation,  A(w,  v;  w;  g)  =  Q(ui)a(w,  v;  n)  + 'H(u>)m(w,v;  n),  avec  Q(u>)  =  1  (resp.,  1  +  iuie),  et 
'H  =  ioj  (resp.,  — w2)  dans  les  cas  parabolique  et  hyperbolique  respectivement ;  g(t)  =  (l/6)f3e_t  est  la 
fonctionnelle  de  controle;  f(v)  est  la  donnee.  L’output  peut  etre  reecrit  comme  dans  (3). 

Nous  introduisons  ensuite  des  espaces  hierarclriques  Xn  de  dimension  N  d’approximation  par  bases 
reduites,  lesquels  sont  generes  par  algorithme  “greedy” .  Une  frequence  ui  et  un  parametre  g,  etant  donnes, 
1’ approximation  par  bases  reduites  satisfait  A(un(u;  g),  v,  w;  g)  =  g(iu>)f(v),\/v  G  Xjj.  L’output  bases 
reduites  est  alors  donne  par  (4).  Enfin,  nous  construisons  un  estimateur  d’erreur  (5).  Comme  enonce  dans 
la  Proposition  1,  la  quantite  A sN(t;n)  borne  l’erreur  entre  l’approximation  par  bases  reduites  (4)  et  la 
solution  elements  finis  de  reference  (3). 

L’output  par  bases  reduites  peut  etre  exprimee  comme  (6).  L’integrale  peut  etre  evaluee  par  residu 
et  est  ecrite  dans  (7).  Dans  cette  formule,  Xn\^)  er  ^n'1  (m)  designent  respectivement  les  fonctions  propres 
et  valeurs  propres  d’un  probleme  aux  valeurs  propres  par  bases  reduites.  Des  procedures  classiques  par 
bases  reduites  offline-online  peuvent  alors  etre  appliquees,  donnant  une  complexite  algorithmique  online 
d’ordre  0(N3  +  Nrif  ).  Des  strategies  similaires  offline-online  peuvent  aussi  etre  employees  pour  le  calcul 
de  la  borne  d’erreur.  La  complexite  algorithmique  online  est  encore  independante  de  la  dimension  A f  de 
l’espace  elements  finis  X-Kf  utilise  pour  le  calcul  de  la  solution  de  reference. 

Les  resultats  numeriques  presentes  dans  la  Figure  1  illustrent  la  precision  et  l’efficacite  de  la  technique 
pour  des  cas  (a)  parabolique,  et  (b)  hyperbolique  :  la  borne  d’erreur  est  petite  (et  l’erreur  reelle  est  encore 
plus  petite),  le  cout  calcul  est  reduit  d’un  ordre  de  grandeur  par  rapport  celui  de  la  solution  de  reference. 


1.  Introduction 

Current  reduced  basis  (RB)  treatment  of  parabolic  equations  [5]  is  quite  effective,  however  RB  treatment 
of  hyperbolic  equations  suffers  from  pessimistic  error  bounds  [12].  Here  we  introduce  a  different  approach 
based  directly  on  continuous  time  rather  than  a  temporal  discretization  —  which  takes  advantage  of 
the  Laplace  transform  (LT)  and  inverse  LT  to  provide  sharper  error  bounds  for  evolution  equations. 

Several  previous  efforts  inform  our  work;  note  here  <r(=  cf>  +  iui)  shall  denote  the  LT  variable.  1) 
Modal  analyses  consider  expansions  in  eigenfunctions  to  yield  reduced  dynamical  systems  [7].  In  our 
work,  we  replace  the  expansion  with  the  inverse  LT,  in  which  a  filter  provides  the  modal  truncation;  we 
replace  the  eigenfunctions  with  the  RB  approximation  of  the  LT  —  for  given  frequency  u>,  a  coercive  or 
noncoercive  elliptic  PDE  [10,13].  Note  we  can  not  provide  rigorous  error  bounds  for  RB  approximations 
of  eigenproblems  [6];  however,  we  can  provide  rigorous  error  bounds  for  the  RB  approximation  of  the 
LT.  2)  The  Krylov/moment-matching  techniques  of  [2,1]  and  Fourier  approaches  of  [3]  consider  the  LT  to 
indentify  a  reduced-order  space;  this  reduced-order  space  then  serves  in  subsequent  Galerkin  projection  in 
the  time  domain.  In  our  work,  we  explicitly  invoke  the  inverse  LT  to  construct  our  RB  approximation  in 
the  time  domain;  this  permits  the  direct  incorporation  of  the  rigorous  RB  (elliptic)  error  bounds  into  the 
parabolic  and  hyperbolic  context.  3)  The  papers  [11,8]  consider  application  of  the  LT/inverse  LT  to  finite 
element  (FE)  semi-discretizations  for  the  purpose  of  parallel  implementation.  In  our  work,  we  accelerate 


2 


the  FE  procedure  by  RB  treatment  of  the  LT;  our  error  bounds  provide  the  necessary  certification.  (We 
also  address  pole-related  issues  through  Online  exact  integration.) 

We  introduce  a  spatial  domain  fl  G  R2  with  boundary  dfl;  we  denote  the  Dirichlet  portion  of  the 
boundary  by  dflD.  We  introduce  the  complex  Hilbert  spaces  L2(fi)  =  {  fn  |c|2  <  oo},  =  {r;  G 

L2(fi)  |  |Vu|  G  L2(fl)},  and  X  =  {v  G  H1( fl)  |u|,9nD  =  0}.  Here  |u|  =  \/vv*  denotes  modulus  and  * 
denotes  complex  conjugate.  We  associate  to  L2(fl)  the  inner  product  (w,v)  =  fnwv*  and  norm  ||ui||  = 
\J (w,  w)  and  to  X  the  inner  product  (w,v)x  =  fn  Vw  •  Vv*  and  induced  norm  ||u?||jc  =  \J (w,  w)x- 
We  introduce  a  real  parameter  /t  which  resides  in  a  closed  bounded  parameter  domain  V  G  Rp.  We 
then  define  parametrized  sesquilinear  forms  m(-,  •;  g)  :  X  x  X  — >  C  (“mass”)  and  a(-,  ■;  g)  :  X  x  X  — >  C 
(“stiffness”);  for  all  /jGD,  m(w,w;g)  and  a(w,w;g)  must  be  real  for  all  w  G  X^ .  We  assume  that  to 
(resp.,  a)  is  symmetric  and  furthermore  continuous  and  coercive  with  respect  to  L2(fl)  (resp.,  X).  We 
also  define  antilinear  bounded  forms  /  :  X  — x  C  (data)  and  £  :  X  — x  C  (output).  Finally,  we  suppose  that 
our  bilinear  forms  to  and  a  are  “affine  in  parameter”  such  that 

Qm  Qa 

m(w,v\n)  =  ^Qqm(g)mq(w,v),  a(w,  v;  g)  =  ^  Qqa(g)aq(w,  v);  (1) 

9=1  9=1 


here  the  O9,  (resp,  0*)  :  V  — x  R,  1  <  q  <  Qm  (resp.,  Qa),  are  /(-dependent  coefficient  functions,  and  the 
mq  (resp.,  a9),  1  <  q  <  Qm  (resp.  Qa),  are  /(-independent  sesquilinear  forms. 

We  now  define  a  “truth”  finite  element  (FE)  approximation  space:  a  standard  P2  polynomial  FE 
approximation  space  X^  C  X  of  dimension  A f .  Our  finite  element  space  X1^  shall  inherit  the  inner 
product  and  norm  associated  to  X.  We  further  define  the  dual  space  (X^)'  and  associated  dual  norm 
II£II(A'a0'  =  sup„eJCjv  |£(u)|/||i>  ||.y-  We  also  introduce  stability  constants  (g)  =  vaiveX^  a(vi  v;g)/ ||u|||- 
and  K^(g)  =  m.iveX*  a(v,  v\  /()/to(i>,  v>  /*)• 

The  parabolic  problem  reads:  Given  any  g  G  V,  find  (the  real  part  of)  u^(t;  g)  G  X ^  such  that 
g),  v;  g)+a(uJ^(t;  g);  v;  g)  =  g(t)f(v),  Vu  G  X^ ,  subject  to  initial  condition  (t  =  0;  /i)  =  0  (for 
simplicity  we  consider  only  homogeneous  initial  conditions) .  We  shall  consider  a  particular  “smooth-start” 
control  function  g(t)  =  (l/6)f3e_t.  The  hyperbolic  problem  reads:  Specify  a  (Rayleigh  or  viscous)  damping 
coefficient,  e  G  R+;  Given  any  g  G  D,  find  (the  real  part  of)  (t;  g)  G  X ^  such  that  m(u£{ (t;  g),v;  g)  + 
ea(u^ (t;  g),v;  g)  +  a(u^(t;  g);  v;  g)  =  g(t)f(v),  \/v  G  ,  subject  to  initial  conditions  w^(t  =  0; /t)  = 
( u^){t  =  0;  g)  =  0.  Our  output  of  interest  (both  for  the  parabolic  and  hyperbolic  cases)  is  then  given  by 


sM(t;g) 


f  B(t  —  t')£(u^ (t';  n))dt' , 
Jo 


(2) 


where  B  is  the  standard  causal  Butterworth  filter  of  order  n f  and  cut-off  frequency  uif.  Note  the  truth 
output  is  defined  by  (2)  and  hence  is  explicitly  filtered:  the  modeler  must  select  rif  and  uj f. 

We  now  state  the  parabolic  and  hyperbolic  problems  in  terms  of  the  LT  and  inverse  LT.  (Note  that  we 
may  consider  here  only  Linear-Time-Invariant  operators/forms.)  We  introduce  a  “combined”  parameter 
g  =  (uj;  g)  which  resides  in  K  x  D  e  Doo-  We  next  define  a  generalized  Helmholtz  problem:  Given 
(uj;g)  G  I?oo5  v^(uj;g)  G  X1^  satisfies  A(u^(uj;  g),  v;  uj;  g)  =  g(iuj)T(v),  \/v  G  X A/”,  where  A(w,v;uj;  g)  = 
G(ui)a(w,  v;  g)  +'H(ui)m(w,  v;  g),  and  X(v)  =  f(v).  (In  the  case  of  non-zero  initial  conditions  there  will  be 
additional  terms  in  F.)  Here  g  is  the  LT  of  g(t):  g(a)  =  1/(ct  +  1)4.  Note  that,  thanks  to  (1),  A(w,  v\  uj;  g) 
admits  an  affine  expansion  in  the  “combined”  parameter  (w;  /t).  We  specify  H  :  R  — X  C  and  Q  :  R  — >  C: 
in  the  parabolic  case,  G{co)  =  1,  'H(w)  =  uj;  in  the  hyperbolic  case,  £/(w)  =  1  +  ioje,  T~L{uj)  =  —uj2. 

We  now  invoke  the  LT  convolution  property  [4]  and  the  inverse  LT  to  express  our  output  (2)  as 


—  /  B(iuj)£{u^ (uj;  g))elultduj. 

2lT  J- 00 


(3) 
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Here  B  is  the  LT  of  the  Butterworth  filter:  for  er  G  C,  B(a)  =  ^fYYk=\{a  ~  af)  \  where  the  a f  = 
u>{  exp(^  +  ^)  exp(  1  <  k  <  rif  ,  are  the  Butterworth  poles. 

In  our  inverse  LT  path  we  have  chosen  the  real  part,  to  be  zero.  In  fact,  a  non-zero  shift  (f)  can 
be  gainfully  exploited.  In  the  parabolic  case,  we  may  choose  negative  <p  (though  still  to  the  right  of  the 
eigenvalues  of  the  differential  operator  and  the  poles  of  the  filter)  to  obtain  decaying  error  bounds.  In  the 
hyperbolic  case,  we  may  choose  a  positive  </>  in  order  to  set  e  =  0  —  no  dissipation;  we  then  obtain  error 
bounds  which  grown  linearly  in  t.  We  treat  these  cases  in  future  work. 


2.  Reduced  Basis  Method 

The  reduced  basis  approximation  shall  be  developed  over  the  “combined”  parameter  /2  =  (w;/Lt);  we 
shall  also  need  the  restricted  parameter  domain  2?  =  [0 ,w]xD  for  some  prescribed  uj  >  oj[  (note  that 
under  our  assumptions  =  (v/^ (ur,  fi))*).  We  first  introduce  the  RB  approximation  spaces  [10] 

relevant  to  both  the  parabolic  and  hyperbolic  case.  We  identify  Nm ax  hierarchical  RB  approximation 
spaces  Xn,  1  <  N  <  iVmax;  here  Xjy  is  of  dimension  N.  These  “Lagrange”  [9,10]  RB  spaces  may  be 
expressed  as  Xn  =  span{£j,i  =  l,...,iV},l  <  N  <  IVmax,  where  the  Q ,  1  <  i  <  Amax,  are  ( i  )x- 
orthonormalized  snapshots  u^(p,Qieed  ),  1  <  i  <  iVmax.  The  sample  points  /XQreed  G  V,  1  <  i  <  iVmax,  at 
which  the  snapshots  are  computed  are  selected  by  a  Greedy  procedure  [13,10].  We  might  also  consider  a 
mixed  approach  with  POD  in  u>  [3]  and  Greedy  in  p  analogous  to  the  time-domain  scheme  of  [5]. 

The  RB  approximation  for  the  LT  is  then  given  by  Galerkin  projection:  Given  (w;/x)  G  Boo,  find 
uat(w;  fi)  G  Xn  such  that  A(un{w;  fi),  v;  w;  p)  =  g(iui)T(v),  V  v  G  Xjq.  Then,  given  t  and  /i,  we  evaluate 
the  RB  output  as 

1  r00  , 

SN(t;n)  =  —  B(iu)£(uN(u-ij,))el“tdu).  (4) 

2n  j- oo 

For  the  appropriate  choice  of  Q  and  T~L  this  formulation  applies  to  both  the  parabolic  and  hyperbolic 
cases.  Galerkin  projection  chooses  a  good  linear  combination  of  the  snapshots. 

We  next  introduce  the  output  error  estimator  for  given  time  t  and  /i  as 

A sN(t;n)  =  27rQ^((AM)^(/t)  (/  (J  \B(iuj)\2\\R(uj; p,)\\2xduj)J  ,  (5) 

where  R  is  the  Riesz  representation  of  the  residual,  (R(lo;h),v)x  =  R(v)  —  giiLu)^1  A{;un{w,  /x),  v;  ui;  /j), 
\/v  G  X-^ ,  and  a is  a  lower  bound  for  cA2  (provided  by  the  Offline-Online  SCM  [10]).  In  the  parabolic 
case  77 (/n)  =  1;  in  the  hyperbolic  case  ?/(/x)  =  r  (e^LB^))1/2)  ’  w^ere  T(z )  —  z(~z/ 2  +  yjl  +  z1  / 4)  and 
is  a  lower  bound  for  KAr(ji)  (constructed  by  variants  on  the  SCM).  We  can  then  state 
Proposition  2.1  For  any  t  >  0  and  fi  GV,  we  obtain  —  sjv(t;/u)|  <  A sN(t-,fi).  □ 

Although  we  only  “train”  the  RB  approximation  over  the  finite  interval  [— uJ,  w],  we  define  the  RB 
approximation  for  all  frequencies;  this  will  be  important  in  order  to  apply  contour  integration.  The 
inaccuracy  of  the  RB  approximation  for  higher  frequency  is  not  of  concern:  the  Butterworth  filter  severely 
attenuates  these  frequencies;  and  we  are  assured  that  the  RB  residual  remains  bounded  thanks  to  stability. 

It  remains  to  develop  a  computational  procedure  for  the  RB  approximation  and  error  bound.  We 
recall  the  Offline-Online  RB  strategy  [10]:  we  admit  significant  Offline  effort  in  exchange  for  greatly 
reduced  cost  in  the  Online  stage  —  in  which  we  aim  to  provide  very  rapid  ( “real-time” )  response  for  each 
new  query  — »  SN{t;ia),  A sN(t;ia).  We  first  introduce  an  RB  eigensystem:  given  p,  a(xiv(/z),  v;  p)  = 

XN(lJ.)m(xN(F),v;iJ.),  \/v  G  XN,  with  eigenpairs  (x^}(/r),  A^}(/Li))  G  (XN,R),  1  <  n  <  N. 

We  may  then  write  the  RB  output  as 


4 


(6) 


N 

sN(t;n)  =  '52i{XN\v))f(x(N)(n))Jn{n) 

n—1 

where  Jn(p)  =  B(iw)g(iw)An(iu;  p)elujtdw;  here  A n{a;  p)  =  1/(A^}(^)  +  CT)  and  A n(cr;/x)  = 

l/((l+ecr)A^(/ti)-|-cr2)  in  the  parabolic  and  hyperbolic  cases,  respectively.  The  integral  is  readily  evaluated 
by  residues  to  yield  (assuming  distinct  poles),  for  the  hyperbolic  problem, 


Jn  =  e»+\ 4n)  -  p(n))-1B(pin))s(pin))  +  /’*(p<n)  -  pi”))-1B(pL”))5(pL”)) 


(7) 


eCTf  *((1  +  e(j}fe)A^  +  a2)  1Bk(&i  )g(<Tf)  +  g  je<T*((l  +  ecr)^A^  +  (J")  1-B(cr)| 


fc  =  l 


where  p±\fi)  =  —eX^\p)/2±  \J —X^\p)  +  (e\^\p) /2)2;  the  parabolic  case  is  similar.  We  observe  the 
connection  to  modal  approaches.  In  the  Online  stage  we  may  assemble  and  solve  the  RB  eigenproblem  in 
0(N3)  operations  and  form  the  f  (Xn) >  £(Xn) >  1  <  n  <  N,  in  0(N2)  FLOPs;  we  then  evaluate  our  output 
(6),  (7)  in  O(Nrif)  FLOPs  per  t  requested.  This  result  relies  on  standard  RB  Offline- Online  procedures 
now  supplemented  with  the  RB  eigenfunction  representation  and  contour  integration. 

The  Offline-Online  approach  for  the  error  bound,  A aN(t;  p)  of  (5),  is  more  involved,  and  the  details 
are  relegated  to  a  future  publication.  The  essential  components  are  the  standard  RB  Offline-Online 
decomposition,  the  RB  eigenfunction  expansion,  and  contour  integration  by  residues.  The  complexity  of 
the  Online  stage  for  the  error  bound  is  an  additional  0(N2(Qm  +  Qa )2  +  N2m). 

We  now  consider  two  model  problems,  one  parabolic  and  one  hyperbolic,  both  posed  on  the  same 
domain  and  over  the  same  P2  FE  truth  approximation  space  of  dimension  J\f  =  9989.  Let  U  be  the  unit 
square  in  R2;  let  Q-|  denote  the  0.5x0. 5  square  centered  at  (0.5,  0.5),  and  define  02  =  The  boundary 

conditions  are  du/dn  =  0  on  the  top  and  bottom  boundaries,  du/dn  =  g(t)  on  the  left  boundary 
and  it  =  0  on  the  right  boundary  dilD .  The  bilinear  forms  are  a(v,w;  p)  =  Vi>  •  Vw  +  p  fn  Vv  ■  Vw, 
m(v,w)  =  JqVW,  and  a(v,w;p )  =  fnVv  ■  Vw,  m(v,w )  =  fQivw  +  p  fnovw,  for  the  parabolic  and 
hyperbolic  problems,  respectively;  the  linear  forms  are  /( v)  =  fgnIttV  and  £{v)  =  f(v );  the  filter  is 
specihed  by  rtf  =  10,  Wf  =  60.  In  the  hyperbolic  case  we  choose  a  damping  coefficient  of  e  =  2E  —  2  (and 
employ  the  lower  bound  Klb(aO  =  (1)/ 1^)-  We  consider  the  parameter  domain  V  =  [1,4]. 

We  generate  a  reduced  basis  with  U  =  120  to  satisfy  an  error  bound  tolerance  etoi  =  IE  —  2  over  a 
Greedy  training  set  of  size  n train  =  10000;  we  require  Amax  =  10  for  the  parabolic  case  and  lVmax  =  85  for 
the  hyperbolic  case.  We  show  in  Figure  1(a)  the  RB  results  for  the  parabolic  case  for  /i  =  1  and  N  =  9; 
the  RB  error  bound  is  A sN(t;  p)  =  6.2 E  —  3  for  all  time  t.  The  Online  RB  (output  and  error  bound)  is  50 
times  faster  than  a  Crank-Nicolson  (At  =  0.1)  FE  truth.  We  show  in  Figure  1(b)  the  RB  results  of  the 
hyperbolic  case  for  p  =  2  and  N  =  75;  the  RB  error  bound  is  A sN(t-,  p)  =  5.7 E  —  3  for  all  time  t.  The 
Online  RB  is  40  times  faster  than  a  second-order  Implicit  Newmark  (At  =  0.05)  FE  truth.  We  introduce 
a  temporal  discretization  of  the  truth  to  more  meaningfully  compare  computational  cost. 
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Figure  1.  Comparison  of  filtered  (cjf  =  60,  rif  =  10)  and  unfiltered  (cjf  =  oo)  FE  solutions  with  the  (filtered)  LT  RB  solution 
for  the  (a)  parabolic  problem,  and  (b)  hyperbolic  problem. 
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